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Abstract —The equations of smoothed particle magnetohy¬ 
drodynamics (SPMHD), even with the various corrections to 
instabilities so far proposed, have been observed to be unstable 
when a very steep density gradient is necessarily combined with 
a variable smoothing length formalism. Here we consider in 
more detail the modifications made to the SPMHD equations in 
LBP2015 that resolve this instability by replacing the smoothing 
length in the induction and anisotropic force equations with an 
average smoothing length term. We then explore the choice of 
average used and compare the effects on a test ‘cylinder-in-a-box’ 
problem and the collapse of a magnetised molecular cloud core. 
We find that, aside from some benign numerical effects at low 
resolutions for the quadratic mean, the formalism is robust as to 
the choice of average but that in complicated models it is essential 
to apply the average to both equations; in particular, all four 
averages considered exhibit similar conservation properties. This 
improved formalism allows for arbitrarily small sink particles 
and field geometries to be explored, vastly expanding the range 
of astronomical problems that can be modeled using SPMHD. 

I. Introduction 

Erom some of the earliest work on SPH by Gingold and 
Monaghan [1] how to include magnetohydrodynamics in the 
equations of SPH has been an area of active research. Whilst 
major progress has been made since then, e.g. the Lagrangian 
formalism derived by Price and Monaghan [2], the correction 
to the tensile instability by Borve et al. [3] (on which there is 
more discussion below), and the divergence cleaning devised 
by [4], SPMHD is known to be unstable in certain astrophys- 
ical simulations. A stable formalism of SPMHD is essential 
to fully study many astrophysical processes, for example the 
generation of protostellar outflows and collimated jets which 
remove angular momentum from an accreting protostellar 
core. 

Extreme gradients in density are not uncommon in astro- 
physical simulations, and these can result in correspondingly 
very small and very large smoothing lengths being needed in 
a physically small region of the simulation. In some cases, 
for example in the pseudo-disc surrounding a protostar the 
smoothing length at the top of one of these gradients can be so 
short that the gradient is not properly sampled, i.e. a situation 
can arise where the particle spacing, Ap In ordinary 

hydrodynamic SPH (and also when radiative transfer schemes 
are employed) no deleterious effects are observed since all the 


equations of SPH either employ both the smoothing length of 
a particle itself and its neighbours or, in the case of the density, 
can be solved self-consistently. However, in the most common 
formalism of magnetohydrodynamics in SPH (see § II) this 
is no longer true, causing a violent instability. We find that 
the instability observed is caused by a very large (and clearly 
unphysical) magnetic field being produced, which then acts 
to accelerate particles rapidly away from the density gradient. 
Importantly, this is not a ‘divergence explosion’ caused by 

a failure to maintain a solenoidal field - observed values of 

/iI-B* I 

' remain <0.1 until after the explosion happens. 

In Lewis, Bate, & Price (submitted) (hereafter ‘LBP2015’) 
we replaced the individual smoothing length terms, 
in these equations with an average term, hab = ^ {ha hb) 
which has the desirable property that at an extreme gradient 
hab tends towards J of the larger value, ensuring the full gra¬ 
dient is sampled and that the equations are evaluated over an 
identical neighbour set. This modified approach is covered in 
more detail in § III. This modified scheme allows for arbitrarily 
small sink particles to be employed, and consequently a much 
larger range of physics to be sampled than hitherto possible. 
We reject simply imposing a minimum smoothing length for 
two reasons, firstly that it is essentially impossible to determine 
a priori a correct value to use; and secondly that in practice 
the minima needed are so large that particle pairing caused by 
stretching the smoothing kernel becomes an serious issue. In 
principle, this could be mitigated by the use of a Wendland 
kernel [5], but is still wasteful of resolution in regions where 
the equations are already stable. The use of an average is more 
nuanced since when ha ^ hb, hab ~ ^{a,6} ^nd therefore 
resolution is not wasted. Additionally, the average can be 
applied to only the equations which are unstable, preserving 
full resolution elsewhere. 

However, the arithmetic mean is not the only plausible 
choice of average - for example the other two Pythagorean 


means both tend to zero if either of h 


{a,6} 


is zero whilst 


the arithmetic mean will approach one-half of the non-zero 
quantity. In § IV we consider the potential advantages of other 
choices of average. In practice, as seen when applied to our test 
‘cylinder-in-a-box’ model in § V, no substantial difference is 
observed. Einally, in § VI we apply these averages to a sample 
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of models similar to those in LBP2015. 


11. Smoothed Particle Magnetohydrodynamics 


We are solving the equations of ideal magnetohydrody¬ 
namics with a gravitational term using the SPMHD method 
initially presented in [2]. We define the MHD stress tensor as 


= -P5^^ + — ( 

Mo V 2 


( 1 ) 


which can be usefully separated into an isotropic ‘pressure’ 
component. 


S^^U = -{P + l—B'^\6^^ 


2 flQ 

and an anisotropic ‘tension’ component, 
5*^'lanis = —BW^ . 


Mo 


We write the total derivative as 




( 2 ) 


( 3 ) 


(4) 


and adopt Einstein’s convention so that repeated indicies imply 
summation. Therefore we can write the equations of ideal 
MHD as 


d 


> 

1 

II 

1 

(5) 

- v> 

dt p 

(6) 

= {B^V^) P - B* , 

(7) 

V‘^4> = AnGp , 

(8) 


where the other symbols have their usual meanings. 

In addition to the discretization in [2], [6], we add artificial 
viscosity and resistivity terms, as a result the final equations 
are not exactly ‘ideal’ MHD. We use the Riemann solver based 
artificial dissipation terms of [7] with spatial and temporally 
varying switches to reduce the dissipation to the minimum 
necessary to maintain numerical stability. The older [8] switch 
is used for the artificial viscosity term, with aAv C [0.1,1.0] 
but the newer [9] switch for artificial resistivity (which is 
observed to exhibit greater stability in protostellar collapse 
simulations) with q^b ^ [0.0,1.0]. The gravitational forces 
are solved using a binary tree and softened using the SPH 
smoothing kernel [10]. 

All magnetic fields found in nature are solenoidal, and 
therefore 

V • B = 0 . (9) 


This constraint is only present in the equations of SPMHD as 
an initial condition; due to round-off error 

( 10 ) 

at 

in general and as a result, the magnetic field will rapidly 
become non-solenoidal and no longer correct. This will inter 


alia produce an unphysical force along the magnetic field lines 
[11] since the term in the anisotropic part of the 

SPMHD momentum equation 

dt p ppo (11) 

=- \{B^V^) B^ + (V^B^)] , 

will no longer be zero. An effective correction to this is to 
subtract a source term exactly equal to this divergence [3], 
which produces an SPMHD momentum equation (neglecting 
the isotropic pressure terms which are unchanged) that de¬ 
pends only on i.e. the smoothing length of each of the 
particles neighbours, viz. 

where Wab is the SPH smoothing kernel and 

O , dhg ^ dWab ( ha ) 

“ dpa V dha 

N 

^ ^ feg A dWab ( K ) 

Vpa ^ dha 


In comparison, the SPMHD induction equation depends 
only upon ha, i.e. the particle’s own smoothing length, 

(14) 

When coupled with a form of magnetic divergence cleaning 
this formalism is remarkably robust. Eor this work, we use 
the constrained hyperbolic divergence cleaning derived in [4], 
where a new scalar field, is coupled to the magnetic field 
such that 

^B%lean = -VV (15) 

and where pj is evolved by 

^ (Y?;') (16) 

where the timescale for damping, 

T= — . (17) 

crCc 

A value of a = 0.8 is used as recommended by [4] to critically 
damp the cleaning wave over a small number of smoothing 
lengths. 


HI. The ‘Average h’ Method 

However, when a very large density gradient is present - e.g. 
in the collapse of a molecular cloud core - this method rapidly 
becomes unstable. If a variable smoothing length regime is 
employed (which is essential in any calculation of this nature) 
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where is a function of p for example, that given in [ 12 ]) 
where 



(where z/ = 3 is the number of spatial dimensions and 77 is 
a parameter controlling the typical number of particles in the 
smoothing sphere) then when Pa ^ ^ ^ ha and 

vice versa it is possible for a particle to have the anisotropic 
component of its magnetic force evaluated over a very large set 
of neighbours and the induction equation evaluated over very 
few. Consequently, particles will interact with each other for 
one equation, but not the other, and an inconsistent estimate 
of both the force and magnetic induction will be calculated. 
Analysis of previous protostellar collapse simulations indicate 
that, for a cubic B-spline kernel with p = 1.2 and therefore 
^ 53 neighbours on average, ratios of greater than 100:4 
neighbours are possible. It is worth noting here that whilst 
77 controls the neighbour count, it does not guarantee that any, 
or every, individual particle will have exactly A^ngh neighbours 
unless the average density profile is flat. 

The result is a violent instability that disrupts the simulation. 
In LBP2015 we replaced the h^a,b} terms in Eqns. (12) 
and (14) with an average term, 

Kb = ^iha + hb) , (19) 

to prevent this instability and were consequently able to 
follow the collapse of a protostar much further than previously 
possible. Whilst this does require the removal of the ^{a,b} 
Eqn. (13) terms, in practice the additional error produced is 
negligible. We therefore obtain in place of Eqn. (12) and 
Eqn. (14), 

d 1 ^ 

-viUs = -^'^{Bi-Bi)BiViWab{hab) , ( 20 ) 
at Ho ^ Pb 

« - 4) BiViWab (Kb). ( 21 ) 

\ Pa / Pa 

IV. Comparison of Averages 


The arithmetic mean used in LBP2015 is not the only 
plausible average. Consequently, we considered the effect of 
choosing a different average. Whilst there are limitless poten- 
tialy viable averages, we only consider the three Pythagorean 
means, viz. the arithmetic mean in Eqn. (19) (hereafter ‘A’), 
the geometric mean (‘G’) defined as 

hab — s/hahi) , ( 22 ) 


and the harmonic mean (‘H’) defined as 


hab 


‘2‘hahi) 
ha h}y 


in addition to the quadratic mean (‘Q’) defined as 



(23) 


(24) 


Table I 

Behaviour of the four means detailed in § IV in limiting cases 


Mean 

ha —>■ 0 

ha = 2hi, 

ha = 10/15 

ha ^ 00 

A 

\hb 

iK 

^hb 

00 

G 

0 

V2hb 

vlo/ 1.6 

00 

H 

0 


fthb 

00 

Q 


Vbhi, 

A-fhb 

00 


All these options take the same value when ha = hi, but 
exhibit significant differences, as seen in table I, in the limiting 
case where ha ^ hi, or vice versa. In particular, whilst A and 
Q are non-zero in the limit h^a,b} ^ 0, G and H are not. It can 
also be shown that for all values of ha and hi,, A > G > H. 
At the other extreme, the differing growth rates of the various 
means can be seen when ha = lO/z^, where Q is nearly 4 
times larger than H. In essence, this affects how much of the 
gradient each average samples, not only in the limit when the 
difference between ha and hi, can be greater a factor of 10 , 
but also when the gradient is shallower. 

V. Numerical Tests 

We use the simple isothermal cylinder-in-a-box test from 
LBP2015 to compare each of these averages and an unmod¬ 
ified code. A cylinder of radius rcyi = 5 code units (units 
defined such that G = 1 and po = 1) with a height-to- 
radius of 2 (2.5 code units thick) and a central hole of radius 
^Tcyi = 0.5 code units was placed in a periodic box with 
a central sink particle to provide a potential equivalent to 10 
times the mass of the material in the cylinder. Sink particles 
(as detailed more comprehensively in [13]) are particles that 
exert no force on the system other than gravity and with an 
‘accretion radius’, r^cc^ whereby any SPH {i.e. gas) particle 
which passes within r^cc of the sink particle is eliminated 
from the simulation and its mass and momentum added to the 
sink particle. We use a sink rather than a simpler potential 
well since this will eliminate any particles which fall out of 
the cylinder and into the centre preventing the timestep from 
becoming needlessly small - since we are using an isothermal 
equation of state, the Courant-limited timestep of particles 
collecting in a central well is very short compared to those in 
the pressure and magnetically supported cylinder. The equation 
of state is given by P (p) = ^up with u fixed so that the 
sound speed was 0.1 code units. An initial magnetic field 
aligned with the z-axis was applied to give a plasma /3, i.e. 
the ratio of hydrodynamic and magnetic pressure, of ^ 8.4. 
(This is equivalent to that derived by assuming the cylinder 
is a sphere of material and using the mass-to-flux equations 
discussed later with p = b). The cylinder was then given a 
differential velocity profile with the initial velocity set to 
obtain a rotation period of T = 2 code units at unit radius. 

We would expect the cylinder material to pile up, forming 
a high density ring with a steep density gradient, so that 
material within a unit radius moves outwards (since it is 
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Figure 1. Density cuts at t = 4.5 of an isothermal test cylinder initially in differential rotation with sink particle providing a central potential equivalent 
to ten times the mass of the cylinder material. The top row (models Dl-4) are the result of using the four averages described in § IV, all of which remain 
stable and conservative. The bottom row shows the unmodified equations (DO) and models where only the induction equation Eqn. (14) (Dla) or anisotropic 
momentum equation Eqn. (12) (Dlb) are modified. A characteristic unphysical bubble can be seen in both DO and Dlb. 


moving faster than the Keplerian velocity) and more distant 
material spiraling inwards. Additionally, some material will 
fall out of the cylinder and towards the sink particle due to 
magnetic and viscous braking effect and the cylinder itself 
will flatten and become more disc like due to rotational and 
self-gravitational forces. 

Calculations were then performed using the four means 
presented above (models Dl-4) and additionally with the 
arithmetic mean but applied to either the induction equation 
only (Dla) or the anisotropic momentum equation only (Dlb). 
We also performed the same simulation with an unmodifled 
code for comparison (DO). 

In Fig. 1 we plot the density profile for all seven models 
just after the explosion has happened in DO. The enhanced 
stability that using an avergae h formalism provides can be 
seen in all four models Dl-4, however, the D4 model (which 
is the quadratic mean) has a somewhat dissimilar profile to 
the three Pythagorean means. The arm-like structures seen in 
some of the plots are caused by numerical artifacts due to 
the very low resolution employed. Models Dla and Dlb show 
that, at least in this simple test, it is sufflcent to apply the 
average to the induction equation only; this is expected since 
the explosion in this instance is clearly driven by an increase in 
magnetic energy as seen in Fig. 2. Even though D4 has evolved 
somewhat differently, no significant difference in conservation 
properties is seen between it and the other three average 
models, as noted in § IV the quadratic mean will produce 
larger values for hab in many cases so the differences seen 
may be due to a somewhat benign loss of effective resolution. 



Eigure 2. The evolution of the total magnetic and kinetic energy for models 
DO (solid black line) and Dl-4 (dashed orange, long-dashed blue, dotted 
green, and solid yellow lines respectively). The timestep shown in Eig. 1 is 
indicated by a vertical line. All four average h models maintain momentum 
conservation whilst the unmodified formalism produces a rapidly growing field 
which ultimately causes an explosion and a consequent increase in kinetic 
energy ca. 1 time unit later. 
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Figure 3. Column density plots at t = 25 670 yrs for simulations of the collapse of a magnetised molecular cloud core of total mass 1 Mq with a 
mass-to-flux ratio of = 5. The top row (models PI-4) are the four average h models described in § IV, the stability and essentially identical evolution 
of which can be clearly seen. Model PO is an unmodified code, which violently explodes. Model Pla is the same as model PI, but with the /i-averaging 
applied only to the induction equation Eqn. (14), whilst this is still stable at this timestep, it becomes unstable shortly after (as seen in the far-right plot on 
the bottom row.) Model Plb has the average applied only to the anisotropic momentum equation Eqn. (12) and, whilst significantly improved over PO, still 
becomes unstable and explodes. 


VI. Protostellar Collapse Models 


The initial conditions are as detailed in LBP2015. In sum¬ 
mary, a IMq sphere of radius 4 x 10^^ cm surrounded by a 
warm medium was placed in a periodic box. The sphere was 
made of ca. 1.5 million SPH particles and the warm medium 
ca. 500,000 particles, which is greater than the requirement to 
resolve the Jeans length in [14]. We set the initial density of the 
sphere to be po = 7.4 x 10“cm“^. The interface between 
the sphere and medium was in pressure equilibrium but with 
a 30:1 density contrast, consequently the initial temperature 
of the sphere is 10 K whilst the surrounding medium is 
approximately 300 K. A barotropic equation state similar to 
that in [15], is used, where 


P = 




P < Pci 

Pci ^ P ^ Pc2 

P > Pc2 


(25) 


to approximate the change in effective 7 (see [16]) as the 
sphere collapses. The two critical densities are set to pd = 
10 “^^g cm“^ ^ lO^po and pc 2 = 10 “^^g cm“^ ^ lO^po- 
We set the initial magnetic field Pq according to the di¬ 
mensionless parameter p, the mass-to-fiux ratio of the sphere. 
This is defined as [17], [18] 



t (yrs) 

Figure 4. The evolution of the total magnetic and kinetic energy for models 
PO (solid black line) and PI-4 (dashed orange, long-dashed blue, dotted green, 
and solid yellow lines respectively). The timestep shown in Fig. 3 is indicated 
by a vertical line. All four average h models approximately maintain energy 
conservation, however, the PO model exhibits a sharp increase in both the 
total magnetic and kinetic energy. Unlike in Fig. 2, the increase in magnetic 
energy does not clearly precede the kinetic energy ‘knee’. 


P = 


Msphere 


/^critical 


(26) 
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where the ratio between the magnetic and gravitational forces 
of the sphere is given by 


/^sphere — 


M 

^’’sphere-So 


(27) 


and the critical ratio (where the magnetic and gravitational 
force are in equilibrium) given by 


_ 2ci r~^ 


(28) 


with Cl =0.53 [19]. We use a mass-to-flux ratio of 5, which 
gives a plasma /3 = 4.5 inside the sphere. The sphere is set in 
solid-body rotation with Q = 1.77 x 10“^^ rad s“^ such that 
the magnitude of the ratio of rotational to gravitational energy 
is initially ^ 0.005. The magnetic field and rotation axes as 
both aligned with the ^-axis, i.e. the parameter defined in 
LBP2015 as the angle between the rotation and field axes, is 
set to = 0°. 

We add sink particles when the density exceeds 
10“^^ g cm“^ ^ 10^ pQ. Previously, accretion radii smaller 
than 5 AU could not be used due to the instability discussed 
earlier (at r^cc > 5 AU the sink is so large that material is 
accreted before a sufficiently large density gradient can be 
created), here we use race = 1 AU but we have also performed 
calculations with significantly smaller sink particles. 

As in § V, we have performed calculations with an unmod- 
ifed code (PO) each average (PI-4), and also with an arithmetic 
mean applied only to the induction equation or momentum 
equation (PIa and Plb). 

Similar to § V, in Fig. 3 we plot the column density of the 
collapsed molecular cloud core att = 25 670 yrs for models 
PO-4. The violent instability in the PO model can be much 
more clearly seen here. We would expect the large cold sphere 
to collapse under its own gravity, and form a thin pseudo¬ 
disc around a dense core, ultimately producing a magnetically 
driven bipolar jet - as seen in PI and in LBP2015. In contrast to 
the earlier tests, the PI a and Plb models both initially appear 
to be stable. However, they eventually both fail indicating 
that in this more complicated model both equations become 
unstable. In Fig. 4 we do not observe the clear separation 
between the increase in magnetic and kinetic energies seen in 
the test model. Unlike the test cylinder model, the quadratic 
mean in P4 does not evolve with significant differences in the 
density profile. 

All four averages produce first hydrostatic core jets with 
velocities of approximately 8 km comparable to those 
obtained with a larger 5 AU sink particle in [20], albeit slightly 
faster due to the smaller sink radius [21]. The evolution of 
these models can then be followed until the jet hits the edge of 
the periodic box ca. 4,000 yrs later. No significant differences 
are seen in the evolution of each of these models, combined 
with the results from models Dl-4 in the previous section 
this indicates that our average h formalism is robust to the 
choice of average provided the full density gradient is sampled, 
however, in some low resolution situations the quadratic mean 
may be a poor choice. 


VH. Conclusion 

The ability to perform stable SPMHD calculations where 
steep density gradients are present has, even with the major 
advances in recent years, been impossible. Having developed 
a slight modification to the equations of SPMHD by using 
an average smoothing length term to fully sample these steep 
density gradients (§ III), we were able in FBP2015 to model 
the collapse of a magnetised cloud core with a < 1 AU 
sink particle and with an arbitrary choice of field geometry. 
Hitherto, such work had been impossible as the calculations 
would violently explode due to (one or more of) a rapidly 
growing field or magnetic force caused by the equations being 
evaluated over extremely dissimilar neighbour sets. In this 
paper we then considered the effect the choice of average 
(§ IV) and which equations are modified on the numerical 
stability of the SPMHD equations. We find that except in low 
resolution tests the formalism is robust to the choice of average 
provided the full density gradient is sampled; however, owing 
to the results in our low resolution tests we recommend the use 
of the Pythagorean means only. Whilst in some test models 
only the SPMHD induction equation needs to be modified to 
ensure stability, in more complicated situations modifications 
to the anisotropic component of the momentum equation is 
also essential to maintain a stable simulation. 
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